Binary Millisecond Pulsar Discovery 
via Gamma-Ray Pulsations 

H. J. Pletsch 1 - 2 *, L. Guillemot 3 , H. Fehrmann 1 ' 2 , B. Allen 1 - 2 ' 4 , M. Kramer 3 ' 5 , C. Aulbert 1 
M. Ackermann 6 , M. Ajello 7 , A. de Angelis 8 , W. B. Atwood 9 , L. Baldini 10 , J. Ballet 11 , 
G. Barbiellini 12 ' 13 , D. Bastieri 14 ' 15 , K. Bechtol 7 , R. Bellazzini 16 , A. W. Borgland 7 , 

E. Bottacini 7 , T. J. Brandt 17 ' 18 , J. Bregeon 16 , M. Brigida 19 ' 20 , P. Bruel 21 , R. Buehler 7 , 
S. Buson 14,15 , G A. Caliandro 22 , R. A. Cameron 7 , P. A. Caraveo 23 , J. M. Casandjian 11 , 

C. Cecchi 24 ' 25 , O. Celik 26 - 27 ' 28 , E. Charles 7 , R.C.G. Chaves 11 , C. C. Cheung 29 , 
J. Chiang 7 , S. Ciprini 30 ' 25 , R. Claus 7 , J. Cohen-Tanugi 31 , J. Conrad 32 ' 33 , S. Cutini 34 , 

F. D'Ammando 24 ' 35 ' 36 , C. D. Dermer 37 , S. W. Digel 7 , P. S. Drell 7 , A. Drlica- Wagner 7 , 

R. Dubois 7 , D. Dumora 38 , C. Favuzzi 19,20 , E. C. Ferrara 26 , A. Franckowiak 7 , 

Y. Fukazawa 39 , P. Fusco 19 ' 20 , F. Gargano 20 , N. Gehrels 26 , S. Germani 24 ' 25 , 
N. Giglietto 19 ' 20 , F. Giordano 19 - 20 , M. Giroletti 40 , G. Godfrey 7 , 1. A. Grenier 11 , 
M.-H. Grondin 41 ' 42 , J. E. Grove 37 , S. Guiriec 26 , D. Hadasch 22 , Y. Hanabata 39 , 
A. K. Harding 26 , P. R. den Hartog 7 , M. Hayashida 7 ' 43 , E. Hays 26 , A. B. Hill 7 ' 44 , 
X. Hou 45 , R. E. Hughes 46 , G. Johannesson 47 , M. S. Jackson 48 ' 33 , T. Jogler 7 , 
A. S. Johnson 7 , W. N. Johnson 37 , J. Kataoka 49 , M. Kerr 7 , J. Knodlseder 17 ' 18 , 
M. Kuss 16 , J. Lande 7 , S. Larsson 32,33,50 , L. Latronico 51 , M. Lemoine-Goumard 38 , 
F. Longo 12,13 , F. Loparco 19,20 , M. N. Lovellette 37 , P. Lubrano 24,25 , F. Massaro 7 , 
M. Mayer 6 , M. N. Mazziotta 20 , J. E. McEnery 26 ' 52 , J. Mehault 31 , P. F. Michelson 7 , 
W. Mitthumsiri 7 , T. Mizuno 53 , M. E. Monzani 7 , A. Morselli 54 , 1. V. Moskalenko 7 , 
S. Murgia 7 , T. Nakamori 49 , R. Nemmen 26 , E. Nuss 31 , M. Ohno 55 , T. Ohsugi 53 , 
N. Omodei 7 , M. Orienti 40 , E. Orlando 7 , F. de Palma 19 ' 20 , D. Paneque 56 ' 7 , 
J. S. Perkins 26 ' 28 ' 27 57 , F. Piron 31 , G. Pivato 15 , T. A. Porter 7 ' 7 , S. Raino 19 ' 20 , 
R. Rando 14 ' 15 , P. S. Ray 37 , M. Razzano 16 ' 9 , A. Reimer 58 ' 7 , O. Reimer 58 ' 7 , T. Reposeur 38 
S. Ritz 9 , R. W. Romani 7 , C. Romoli 15 , DA. Sanchez 41 , P. M. Saz Parkinson 9 , 
A. Schulz 6 , C. Sgro 16 , E. do Couto e Silva 7 , E. J. Siskind 59 , D. A. Smith 38 , 

G. Spandre 16 , P. Spinelli 19 ' 20 , D. J. Suson 60 , H. Takahashi 39 , T. Tanaka 7 , J. B. Thayer 7 , 

J. G. Thayer 7 , D. J. Thompson 26 , L. Tibaldo 14 - 15 , M. Tinivella 16 , E. Troja 26 , 
T L. Usher 7 , J. Vandenbroucke 7 , V. Vasileiou 31 , G. Vianello 7 ' 61 , V. Vitale 54 ' 62 , 
A. P. Waite 7 , B. L. Winer 46 , K. S. Wood 37 , M. Wood 7 , Z. Yang 32 ' 33 , S. Zimmer 32 ' 33 

*To whom correspondence should be addressed. Email: holger.pletsch@aei.mpg.de 
Affiliations are listed at the end of the paper. 



1 



Millisecond pulsars, old neutron stars spun-up by accreting matter from a companion star, can 
reach high rotation rates of hundreds of revolutions per second. Until now, all such "recy- 
cled" rotation-powered pulsars have been detected by their spin-modulated radio emission. In a 
computing-intensive blind search of gamma-ray data from the Fermi Large Area Telescope (with 
partial constraints from optical data), we detected a 2.5-millisecond pulsar, PSR J1311 3430. 
This unambiguously explains a formerly unidentified gamma-ray source that had been a decade- 
long enigma, confirming previous conjectures. The pulsar is in a circular orbit with an orbital 
period of only 93 minutes, the shortest of any spin-powered pulsar binary ever found. 



Almost exactly 30 years ago, radio observations de- 
tected the first neutron star with a millisecond spin pe- 
riod (7). Launched in 2008, the Large Area Telescope 
(LAT) on the Fermi Gamma-ray Space Telescope (2) 
confirmed that many radio-detected millisecond pulsars 
(MSPs) are also bright gamma-ray emitters (3). In each 
case, gamma-ray (0.1 to 100 GeV) pulsations were re- 
vealed by using rotation parameters obtained from ra- 
dio telescopes (4) to assign rotational phases to LAT- 
detected photons. 

The Fermi LAT also provides sufficient sensitivity 
to detect pulsars via direct searches for periodicity in 
the sparse gamma-ray photons. Such blind searches (5) 
of LAT data for solitary pulsars have so far unveiled 
36 younger gamma-ray pulsars (6-9) with rotation rates 
between 2 and 20 Hz. In the radio band, all but four 
of these objects remain completely undetected despite 
deep follow-up radio searches (10). This is a large frac- 
tion of all young gamma-ray emitting neutron stars and 
shows that such blind-search gamma-ray detections are 
essential for understanding the pulsar population (11). 
However, no MSP has been detected via gamma-ray 
pulsations until now, and so we have not been able to 
see whether a similar population of radio-quiet MSPs 
exists. 

The blind-search problem for gamma-ray pulsars is 
computationally demanding, because the relevant pul- 
sar parameters are unknown a priori and must be ex- 
plicitly searched. For observation times spanning sev- 
eral years, this requires a dense grid to cover the multi- 
dimensional parameter space, with a tremendous num- 
ber of points to be individually tested. Blind searches 
for MSPs in gamma-ray data are vastly more difficult 
than for slower pulsars largely because the search must 
extend to much higher spin frequencies [to and beyond 
716 Hz (12)]. Furthermore, most MSPs are in binary 
systems, where the additionally unknown orbital pa- 
rameters can increase the computational complexity by 
orders of magnitude. Thus, blind searches for binary 
MSPs were hitherto virtually unfeasible. 



We have now broken this impasse, detecting a bi- 
nary MSP, denoted PSR J131 1-3430, in a direct blind 
search of the formerly unidentified gamma-ray source 
2FGL J1311.7— 3429, one of the brightest listed in the 
Fermi-LAT Second Source Catalog [2FGL (73)]. This 
source also had counterparts in several earlier gamma- 
ray catalogs and was first registered in data from the 
Energetic Gamma Ray Experiment Telescope [EGRET 
(14)] on the Compton Gamma Ray Observatory. 

In a search for potential optical counterparts of 
2FGL J1311.7— 3429, Romani (75) identified a quasi- 
sinusoidally modulated optical flux with a period of 
93 minutes and conjectured this to be a "black-widow" 
pulsar binary system (16). In this interpretation, an MSP 
strongly irradiates what is left of the donor companion 
star to eventually evaporate it. This plausibly explained 
the observed brightness variation resulting from strong 
heating of one side of the companion by the pulsar ra- 
diation. Associating this optical variation with the or- 
bital period of the putative binary system constrained 
the ranges of orbital search parameters and also con- 
fined the sky location for the search. Thus, these con- 
straints made a blind binary-MSP search in LAT data 
feasible; however the computational challenge involved 
remained enormous. To test the binary-MSP hypothesis 
as the possible nature of 2FGL J1311.7— 3429, we de- 
veloped a method to search the LAT data for pulsations 
over the entire relevant parameter space. 

Under the black-widow interpretation, the search is 
confined toward the sky location of the potential opti- 
cal counterpart and the orbit is expected to be circu- 
lar, leaving a five-dimensional search space. The in- 
dividual dimensions are spin frequency /, its rate of 
change /, the orbital period P or b, time of ascending 
node T asc , and x — a p sin i, the projection of the pul- 
sar semi-major axis a p onto the line of sight with or- 
bital inclination angle i. We designed the blind search 
to maintain sensitivity to very high pulsar spin frequen- 
cies, / < 1.4 kHz, and values of / typical for MSPs, 
-5 x 10 _13 Hzs _1 < / < 0. Although the optical 
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data constrain P or b and T asc , the uncertainties are by 
far larger than the precision necessary for a pulsar de- 
tection. This required us to search ranges of P or b = 
5626.0 ± 0.1s and T asc = 56009.131 ± 0.012 MJD 
(modified Julian days) around the nominal values (75), 
and < x < 0.1 lt-s (light-seconds). 

Searching this five-dimensional parameter space 
fully coherently given a multiple-year data time span 
is computationally impossible. To solve this problem, 
we used the hierarchical (three-staged) search strategy 
that previously enabled the detection of 10 solitary, 
younger (i.e. non-MSP) pulsars in blind searches of 
LAT data (8, 9), exploiting methods originally devel- 
oped to detect gravitational waves from pulsars (17-20). 
Here we expanded this approach to also search over bi- 
nary orbital parameters. The first stage of the hierar- 
chical scheme is the most computing-intensive and uses 
an efficient "semi-coherent" method (8), extending the 
method of Atwood et al. (27). This step involves (inco- 
herently) combining coherent Fourier power computed 
using a window of 2 20 s (~12 days) by sliding the win- 
dow over the entire LAT data set (hence the term "semi- 
coherent"). In a second stage, significant semi-coherent 
candidates are automatically followed up through a fully 
coherent analysis made possible because only a small 
region of parameter space around the candidate is ex- 
plored. A third stage further refines coherent pulsar can- 
didates by including higher signal harmonics [using the 
77-test (22, 23)]. The computing cost to coherently fol- 
low up a single semi-coherent candidate is negligible 
relative to the total cost of the first stage. Therefore, 
constructing the search grid of the semi-coherent stage 
as efficiently as possible is of utmost importance. 

The key element in constructing an optimally effi- 
cient grid for the semi-coherent search is a distance met- 
ric on the search space (17-19,24). The metric provides 
an analytic geometric tool measuring the expected frac- 
tional loss in signal-to-noise ratio (S/N) squared for any 
given pulsar-signal location at a nearby grid point. The 
metric is obtained from a Taylor expansion of the frac- 
tional loss to second order around the parameter-space 
location of a given signal. In contrast to searching for 
solitary pulsars, a difficulty in the binary case is that the 
metric components explicitly depend on the search pa- 
rameters (24). Thus, the metric (and so the grid-point 
density required to not miss a signal) changes across 
orbital parameter space. Constructing a simple lattice 
with constant spacings would be highly inefficient, re- 
sulting in either vast over- or under-covering of large 



parameter-space regions. We developed a grid construc- 
tion algorithm (25) that effectively uses the metric for- 
malism. Orbital grid points were first placed at random, 
then those that were either too close together or too far 
apart according to the metric were moved (barycentric 
shifts), minimizing the maximum possible loss in S/N 
for any pulsar signal across the entire search parameter 
space. By design the resulting grid (25) ensured never 
losing more than 30% in S/N for any signal parameters. 

The input LAT data we prepared for this search 
spanned almost 4 years (1437 days) and includes 
gamma-ray photons with LAT-reconstructed directions 
within 15° around the targeted sky position (25). To im- 
prove the S/N of a putative pulsar signal, we assigned 
each photon a weight (23) measuring the probability of 
originating from the conjectured pulsar, computed with 
a spectral likelihood method (25). The gamma-ray spec- 
trum of 2FGL J1311.7— 3429 is best modeled by an ex- 
ponentially cut-off power law (Fig. SI), with spectral 
parameters reminiscent of other gamma-ray pulsars (Ta- 
ble 1). The computational work of the search was done 
on the ATLAS cluster in Hannover, Germany. Soon af- 
ter initiation, the searching procedure convincingly de- 
tected PSR J 131 1-3430. 

Following the blind-search detection, we refined the 
pulsar parameters further in a timing analysis (26). We 
obtained pulse times of arrival (TOAs) from subdi- 
viding the LAT data into 40 segments of about equal 
length. We produced a pulse profile for each segment 
using the initial pulsar parameters, and cross-correlated 
each pulse profile with a multi-Gaussian template de- 
rived from fitting the entire data set to determine the 
TOAs. We used the Tempo2 software (27) to fit the 
TOAs to a timing model including sky position, /, /, 
and binary-orbit parameters (Fig. 1 and Table 1). We 
found no statistically significant evidence for orbital 
eccentricity at the e < 10~ 3 level. We measured a 
marginal evidence for a total proper motion of 8 ± 3 
milliarcseconds/year. Generally, the observed value of 
/ = (-3.198 ±0.002) x lO^Hzs- 1 is only an upper 
limit of the intrinsic frequency change /i n , because of 
the Shklovskii effect in which Doppler shifts caused by 
the proper motion can account for part of /. Under the 
assumption that the proper motion of PSR J 13 1 1 —3430 
is small enough to approximate / = f 1D , we derived fur- 
ther quantities from the pulsar rotational parameters (Ta- 
ble 1). 

The rotational ephemeris of PSR J131 1—3430 also 
provides constraints on the companion mass m c through 
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the binary mass function that combines x, P m b, and the 
gravitational constant G, 

47T 2 x 3 (m c sin/,) 3 



= (2.995 ± 0.003) x l(r 7 M (1) 

where m p is the pulsar mass and M Q is the mass of 
the Sun. Typical MSP masses are 1.35 to 2.0 M . 
Assuming m p = 1.35 M and i = 90° (orbit is 
edge-on) yields the minimum companion mass, m c > 
8.2 x 10~ 3 M Q , which is only about eight times the 
mass of Jupiter. By means of Kepler's third law and 
typical MSP masses (m p >• m c ), the binary sepa- 
ration, a = a p + a c is accurately approximated by 
a = 0.75R Q (mp/1.35MQ) 1/3 , where R is the radius 
of the Sun. Thus, PSR J131 1-3430 is likely the most 
compact pulsar binary known. 

The compact orbit and the optical flaring events (75) 
suggest that the pulsar heating is driving a strong, possi- 
bly variable, stellar wind of ablated material of the com- 
panion. Interactions with the companion wind could af- 
fect the gamma-ray flux observed. In a dedicated anal- 
ysis (25), we found no evidence for a modulation at the 
orbital period of the gamma-ray flux or its spectrum. 

We also examined the gamma-ray spectral param- 
eters of PSR J131 1—3430 as a function of rotational 
phase (25). Dividing the data into 10 segments ac- 
cording to different rotational-phase intervals, we spec- 
trally analyzed each segment separately. In line with the 
background estimation in the pulse profile (Fig. 1), we 
detected significant gamma-ray emission at all phases. 
The gamma-ray spectrum in the off-pulse phase inter- 
val (Fig. S2) is better modeled by an exponentially cut- 
off power law, potentially indicative of magnetospheric 
origin from the pulsar, rather than by a simple power 
law which would more likely suggest intrabinary wind 
shock emission. 

Repeated, sensitive radio searches of the previously 
unidentified gamma-ray source, including Green Bank 
Telescope observations at 820 MHz gave no pulsar de- 
tection (28). However, material ablated from the com- 
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panion by the pulsar irradiation might obscure radio 
pulses. At higher radio frequencies decreased scattering 
and absorption resulting in shorter eclipses are observed 
for other black-widow pulsars (72). 

The optical observations provide evidence for strong 
heating of the pulsar companion that is near filling its 
Roche lobe (75). With m p = 1.35 M and i = 90°, 
the Roche lobe radius of the companion is to good ap- 
proximation (29) Rl = 0.063 Rq. The minimum mean 
density of the Roche-lobe filling companion directly fol- 
lows from the orbital period (30), p = 45gcm~ 3 . This 
is twice the density of the planetary-mass companion of 
PSR J1719-1438 (37). One scenario for the formation 
of that system posits an ultra-compact x-ray binary with 
a He or C degenerate donor transferring mass to the neu- 
tron star. However, van Haaften et al. (32) argue that 
angular momentum losses through gravitational-wave 
emission are insufficient to reach the low masses and 
short period of the PSR J1719— 1438 system within the 
age of the universe. Instead, strong heating to bloat 
the companion or extra angular momentum loss from 
a companion evaporative wind are required. An alterna- 
tive scenario (33) proposes that a combination of angu- 
lar momentum loss and wind evaporation from an initial 
companion mass of 2M in a 0.8 day orbit can bring 
the system to low masses and short orbital periods in 
~ 6 billion years. Indeed, their scenario produces a 
good match to the m c <~ 0.01M Q , P orb ~ 0.065 days 
seen for PSR J131 1-3430. At this point in the evolution 
the system is detached, the companion is He-dominated 
and irradiation has taken over the evolution. Presum- 
ably continued irradiation can drive the system towards 
PSR J 1719 1438-type companion masses, or produce 
an isolated MSP. 

The direct detection of an MSP in a blind search 
of gamma-ray data implies that further MSPs, includ- 
ing other extreme binary pulsars, may exist among the 
bright, as-yet unidentified 2FGL gamma-ray sources 
[e.g. (34, 35)] which are too radio faint or obscured by 
dense companion winds to be found in typical radio 
searches. 
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Table 1. Measured and derived parameters for PSR J 1 3 1 1 — 3430, with formal la uncertainties in the last digits 
(dd, days; hh, hours; mm, minutes; ss, seconds). Spectral parameters are averages over pulse phase. 



Parameter 


Value 


Right ascension (J2000.0) (hh:mm:ss) 


13:11:45.7242(2) 


Declination (J2000.0) (dd:mm:ss) 


-34:30:30.350(4) 


Spin frequency, / (Hz) 


390.56839326407(4) 


Frequency derivative, / (10~ 15 Hz s" 1 ) 


-3.198(2) 


Reference time scale 


Barycentric Dynamical Time (TDB) 


Reference time (MJD) 


55266.90789575858 


Orbital period P or b (d) 


0.0651157335(7) 


Projected pulsar semi-major axis x (lt-s) 


0.010581(4) 


Time of ascending node T asc (MJD) 


56009.129454(7) 


Eccentricity e 


< 0.001 


Data span (MJD) 


54682 to 56119 


Weighted RMS residual (/xs) 


17 


Derived Quantities 


Companion mass m c (M Q ) 


> 0.0082 


Spin-down luminosity E (erg s" 1 ) 


4.9 x 10 34 


Characteristic age r c (years) 


1.9 x 10 9 


Surface magnetic field B§ (G) 


2.3 x 10 8 


Gamma-Ray Spectral Parameters 


Photon index, T 


1.8(1) 


Cutoff energy, £ c (GeV) 


3.2(4) 


Photon flux above 0.1 GeV, F (10~ 8 photons cm" 2 s" 1 ) 


9.2(5) 


Energy flux above 0.1 GeV, G (10" n erg cm" 2 s" 1 ) 


6.2(2) 
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Fig. 1. Phase-time diagram and gamma-ray pulse profiles for PSR J 131 1—3430. Two pulsar rotations are shown 
for clarity. (Left) The pulsar rotational phase for each gamma-ray-photon arrival time; probability weights are 
shown in color code. (Right) The pulse profiles in different energy bands. Each bin is 0.02 in phase and photon 
weights are used. The dashed line indicates the estimated background level from a surrounding annular region. 
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Supporting Online Material 



Materials and Methods: Fermi-LAT Data Analysis 



LAT data selection 

The LAT surveys the entire sky every 3 hours (two or- 
bits). In this work, we used data taken in this sky-survey 
mode between 4 August 2008 and 10 July 2012. The 
data were processed using the Fermi Science Tools 1 
(v9r28p0). We selected gamma-ray photons belonging 
to the "Source" class under the P7V6 event selections. 
Photons with reconstructed zenith angles larger than 
100° were rejected in order to exclude the bright contri- 
bution from the Earth's limb. We also rejected photons 
recorded when the instrument was not operating in sky- 
survey mode, or when its rocking angle exceeded 52°. 
Finally, we selected only photons with energies > 0.1 
GeV and found within 15° of the direction of the pulsar. 

LAT data spectral analysis 

We determined the gamma-ray spectral properties of 
PSR J131 1—3430 by performing a binned likelihood 
analysis of this data set, using the pyLikelihood 
tool. The Galactic diffuse emission was modeled us- 
ing the glLiem_v02 map cube, and the extragalactic 
emission and the residual instrumental backgrounds 
were modeled jointly using the iso_p7v6source tem- 
plate 2 . The spectral model used in this analysis also 
included the contributions of 2FGL sources (13) within 
20° of the center of the field of view, and the contri- 
bution from the pulsar was modeled using an exponen- 
tially cut-off power-law (ECPL) of the form dN/dE oc 
E~ T exp(— E/E c ) where E represents the photon en- 
ergy, T is the photon index and E c is the cutoff en- 
ergy of the spectrum. The normalizations of the dif- 
fuse components were left free in the fit, as well as 
spectral parameters of sources within 8° of the pul- 
sar sky position. The best-fit spectral parameters for 
PSR J13 11-3430 along with the derived photon and 
energy fluxes are displayed in Table 1 and character- 
ize the pulse-phase-averaged gamma-ray spectrum of 
the pulsar (Fig. SI). The measured energy flux is con- 
sistent with that published earlier for this gamma-ray 
source (73). We checked the best-fit parameters listed 
in Table 1 with the point like likelihood analysis 
tool (23), and found results that are consistent within 

'http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html 
2 http://fermi.gsfc.nasa.gOv/ssc/data/access/lat/BackgroundModels.l 



uncertainties. The tool gtsrcprob and the spectral 
model obtained from the likelihood analysis were finally 
used to calculate the probabilities that the photons orig- 
inate from PSR J131 1-3430. 

We also performed a spectral analysis re- 
solved in pulse phase (rotational phase) for 
PSR J131 1-3430. For this, we divided the data set into 
ten segments according to pulse phase and measured 
the gamma-ray spectrum in each of those segments in 
a binned likelihood analysis assuming an ECPL model 
for the pulsar (Fig. S2). As observed for other bright 
gamma-ray pulsars [e.g. (36)], the spectral properties 
of PSR J13 11-3430 evolve strongly with rotational 
phase, suggesting varying emission altitudes and cur- 
vature radii of the magnetic field lines. In addition, 
significant gamma-ray emission is detected over all ro- 
tational phases: selecting phases between 0.8 and 1.2 
we measured a photon flux for PSR J 131 1—3430 of 
(4.6 ±0.9) x 10~ 8 photons cm -2 s~\ and a photon flux 
of (2.8±1.6) x 10~ 8 photons cm~ 2 s _1 forpulse phases 
between 0.9 and 1.1. In the former phase interval, the 
exponentially cut-off power-law model is preferred over 
a simple power-law spectral shape at the <~ 3.5cr signif- 
icance level, and at <~ 2.5a for the latter phase interval. 
This potentially indicates a magnetospheric origin for 
this "off-pulse" emission. The existence of magneto- 
spheric emission in the off-pulse region of the gamma- 
ray pulse profiles is predicted from theoretical models 
under specific geometrical orientations (37) and can 
give insights into the pulsar emission geometry. 

We also conducted an analysis to look for an orbital 
modulation of the gamma-ray flux and the spectrum. We 
subdivided the orbit into ten equally spaced bins, and for 
each subset of data we performed a binned likelihood 
analysis. The flux, the spectral index T, and the energy 
cutoff E c are compatible with a constant value all along 
the orbit. Thus, we proceeded to compute the formal 
95% confidence upper limits on the amplitude of a sinu- 
soidal modulation and obtained < 2.6 x 10~ 8 photons 
cm~ 2 s" 1 for the flux, < 0.32 for T, and < 2.2 GeV 
for E c . 



11 



Pulsar search in LAT data 

To correct for the Doppler modulation due to the satel- 
lite's motion in the Solar System, we applied barycenter 
corrections to the arrival times of the LAT gamma-ray 
photons using the JPL DE405 Solar System ephemeris. 

Constructing the parameter-space grid for the semi- 
coherent search (the first stage of the hierarchical search 
scheme) as efficiently as possible is of utmost impor- 
tance, because this stage dominates the overall comput- 
ing cost. For this purpose, we developed an algorithm 
that effectively utilizes the metric formalism. The met- 
ric provides a geometric tool measuring the expected 
fractional loss in squared signal-to-noise ratio for any 
given pulsar-signal location at a nearby grid point. The 
metric components along the search-space directions of 
/ and / are constant across the entire space [see, e.g., 
(18, 19)]. In contrast, the orbital metric components (in 
search-space directions of {P or b, T asc , x}) explicitly de- 
pend upon the search parameters [see, e.g., (24)]. This 
implies that the metric (and so the grid-point density re- 
quired to not miss a signal) changes across orbital pa- 
rameter space. Therefore, constructing a simple lattice 
with constant spacings over these dimensions would be 
highly inefficient, resulting in either vast over- or under- 
covering of large parameter-space regions. 

In contrast, the grid-construction algorithm we de- 
veloped follows a more efficient approach. While us- 
ing constant spacings in / and /, grid points over 
{P rb,T asc ,x} are first placed at random [e.g., (38)]. 
Then those that are either too close together or too 
far apart according to the metric are moved (barycen- 
tric shifts), minimizing the maximum possible loss in 
signal-to-noise ratio for any pulsar signal across the en- 
tire search parameter space. To accelerate this otherwise 
computationally-bound process, we divided the search 
space into sub-volumes to avoid metric distance com- 
parisons between a trial point and every other grid point, 
exploiting an efficient hashing technique. For an arbi- 
trary point in parameter space, the index of the enclos- 
ing sub-volume (and thus the parameters of its neigh- 
boring points) is obtained from a hash table through a 
fast rounding operation. By design the resulting grid en- 
sured never losing more than 30% of the signal-to-noise 
ratio for any signal parameters. Finally, we used sim- 
ulated pulsar signals to validate this design goal. This 
process is also highly accelerated by exploiting the hash 
table, because it provides quick access to the closest 



grid points around any given parameter-space location. 
This way, the search for each simulated pulsar signal 
is computationally inexpensive, because only the rele- 
vant subset of nearest grid points around each signal are 
searched. 

In the blind search, we obtained the semi-coherent 
detection statistic (coherent Fourier power) computed 
using a coherence window of T — 2 20 s (~12 days) be- 
ing incoherently combined by sliding the window over 
the entire 4 years of LAT data) over the entire / grid 
by exploiting the efficiency of the fast Fourier transform 
(FFT) algorithm. For this purpose, we divided total / 
search range into separate bands using a heterodyning 
bandwidth of A/bw = 128 Hz. Thus, the FFT contains 
T A/bw — 10 8 frequency bins. This choice for A/bw 
allowed us to fit the computation into memory on the 
ATLAS computing cluster to maximize performance. 
In the / direction, we analyzed about 10 2 uniformly 
spaced grid points. The use of separate frequency bands 
is also extremely advantageous in view of the orbital 
grids, which we adapted to each band. This further re- 
duced the computational cost, since the total number 
of required grid points to cover the orbital parameter 
space {P rb, Paso x} is about 10 7 (/ max /700Hz) 3 and 
increases with / max cubed, where / max is the highest 
spin frequency (most conservative choice) in the search 
band (before heterodyning). Thus, in total the search 
grid covering the frequency band near 700 Hz com- 
prised about 10 8 x 10 2 x 10 7 = 10 17 points. We note that 
the heterodyning step reduces the bandwidth of the data 
to 128 Hz, which is still much larger than the narrow 
bandwidth of the signal (that is about / x 10~ 4 , which 
is - 0.04 Hz for PSR J131 1-3430) that remains unaf- 
fected by this. The heterodyning does not alter the pul- 
sar signal shape (i.e., the extension in parameter space) 
- it is merely shifted in the frequency domain and so the 
metric is unaffected and still depends on / max , the orig- 
inal upper frequency. This /^ ax -dependency is read- 
ily seen from the fact that the amplitude of the pulsar- 
rotational-phase modulation due to the binary motion is 
proportional to the spin frequency /. And, the met- 
ric components involve products of first-order deriva- 
tives of the rotational phase with respect to the param- 
eters (17-19, 24). As the number of grid points is pro- 
portional to the square root of the metric determinant, 
it therefore increases with / 3 (one / contribution from 
each of the three orbital dimensions searched). 
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Fig.Sl. Pulse-phase-averaged gamma-ray spectral energy distribution for PSR J 1 3 1 1 — 3430. Data points (solid 
squares) are derived from fits of individual energy bands with variable widths, in which the pulsar is detected with 
a significance greater than 15a. In these individual bands the pulsar emission is modeled with a simple power-law 
spectrum. We calculated an upper limit for the highest energy band, in which the pulsar was not detected with 
sufficient significance. The solid curves represents the best-fit model obtained from the likelihood analysis and 
dashed curves indicate the la uncertainties. 
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Fig.S2. Pulse-phase-resolved gamma-ray spectral analysis for PSR J131 1—3430 using 10 bins per rotation. The 
different panels show the photon index (top left), cutoff energy (top right), photon flux above 0. 1 GeV (bottom 
left), and energy flux above 0.1 GeV (bottom right). Error bars indicate statistical la uncertainties. 
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